Using Householder Matrices to 
Establish Mixing Test Critical Values 



Abstract: A measure-preserving dynamical system can be approxi- 
mated by a Markov shift with a bistochastic matrix. This leads to using 
empirical stochastic matrices to measure and estimate properties of stir- 
ring protocols. Specifically, the second largest eigenvalue can be used 
to statistically decide if a stirring protocol is weak-mixing, ergodic, or 
nonergodic. Such hypothesis tests require appropriate probability distri- 
butions. In this paper, we propose using Monte Carlo empirical probabil- 
ity distributions from unistochastic matrices to establish critical values. 
These unistochastic matrices arise from randomly constructed House- 
holder matrices. 
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1. Introduction 

If a dynamical system has a probability measure and the system is measure- 
preserving, then partitioning the domain into n states of equal measure leads 
to a Markov shift whose measure is defined by a bistochastic matrix and the 
length n row vector 

(^••^)- (1-1) 

This partition approximation is called Ulam's method [21]. After partitioning 
the space, data point movement from one iteration of the function provides 
a stochastic matrix that approximates the bistochastic matrix from Ulam's 
method. Hence a Markov shift with an empirical stochastic matrix and the 
1/n row vector approximates the dynamical system [see [17, Chapter 1, Chap- 
ter 9] for procedure and convergence rate]. 

We may model a stirring protocol's affect on a compression-resistant fluid 
with a measure-preserving dynamical system. In this paper we are interested 
in discrete interations of a stirring protocol where the fluid at the beginning 
is the same fluid at the end. We 'look' at the fluid before and after stirring, 
but not during. 

Properties of an empirical Markov shift can measure and evaluate proper- 
ties of a measure-preserving dynamical system [8, 10, 13]. The second largest 
eigenvalue of a empirical stochastic matrix arising from Ulam's method may 
be used to statistically decide if a measure-preserving dynamical system is 
weak-mixing, ergodic, or nonergodic [6, 7, 9, 14, 17]. To statistically test if 
the dynamical system is ergodic, we need to have some knowledge of 

P(|A2-1|>A;:A2 = 1); (1.2) 

to statistically test if the dynamical system is weak-mixing, we need to have 
some knowledge of 

PdAsI > : IA2I = 1). (1.3) 

We use A2 to denote the second largest eigenvalue of a bistochastic matrix 
arising from Ulam's method, and A2 denotes the second largest eigenvalue 
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of a corresponding empirical stochastic matrix (bistochastic and unistochas- 
tic matrices will be defined shortly). The utility of A2 as a test statistic 
arises directly from the relationship between stochastic matrix eigenvalues 
and Markov shift ergodic, mixing properties. Since the unit circle contains 
all eigenvalues of stochastic matrices, there are no reasonable probability 
distributions of A2 with 1 as the mean or median of either A2 or IA2I. So for 
hypothesis testing, we should use a probability distribution that has signifi- 
cant mass near A2 = 1 or IA2I = 1. 

In this paper, we show that it is reasonable to approximate the condi- 
tional probability distributions with Monte Carlo probability distributions 
when the equal-measure partition sets are small. These Monte Carlo prob- 
ability distributions are constructed using randomly generated Householder 
matrices. 

Stirring protocols of compression resistant fluids, such as chocolate and 
water, provide examples of nearly measure-perserving dynamical systems. 
The need for confidence in the mixing of food items and in the mixing of 
pharmaceuticals highlights the utility of such probability distributions. 

We propose using randomly generated, nonzero, independent, indentically 
distributed real numbers to generate Householder matrices; take products of 
permutation matrices with Householder matrices; then square the magnitude 
of the products' entries to get unistochastic matrices. From these unistochas- 
tic matrices, construct a Monte Carlo approximation of a desired probability 
distribution. From the Monte Carlo probability distribution, establish the 
critical value for rejecting the null-hypothesis. The primary focus of this pa- 
per is to establish a method for determining hypothesis test critical values. 
Deciding which specific probability distribution to use in a hypothesis test 
depends on properties of the dynamical system; we will only show that the 
presented Monte Carlo methods are reasonable and leave probability distri- 
bution selection for the future. 

There are several ways to use Monte Carlo methods to generate bistochas- 
tic matrices; unfortunately, many techniques lead to empirical probability 
distributions where the central tendency of A2 is close to zero [17, Chapter 
12]. Such distributions provide little utility for a weak-mixing, ergodic, or 
nonergodic hypothesis test. The methods presented here lead to Householder 
matrices that, in a Frobenius norm sense, are likely to be close to the identity 
matrix. Squaring the magnitude of entries from these unitary matrices gives 
unistochatic matrices. If we want a unistochastic matrix close to a particular 
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permutation matrix, we may multiply the Householder matrix by the desired 
permutation matrix. The main advantage of this method is that it provides 
probability distributions based on observed unistochastic matrices. 

2. Bistochastic Matrices 

Definition 2.1. An n x n bistochastic matrix is a stochastic matrix whose 
transpose is also a stochastic matrix. 

By the Birkhoff-von Neumann theorem, the set of n x n bistochastic ma- 
trices form a convex set with permutation matrices as extreme points. We 
refer to this set as Birkhoff's polytope [2, 3]. Bistochastic matrices are also 
referred to as doubly stochastic. 

Definition 2.2. An n x n bistochastic matrix is called unistochastic if each 
entry is equal to the squared magnitude of some unitary matrix. 

The set of n x n unistochastic matrices form a proper subset of Birkhoff's 
polytope [2, page 307, section 1]. Since the set of unistochastic matrices is a 
proper subset, the proposed method should only be used when the Ulam's 
method bistochastic matrix is approximately unistochastic. 

Definition 2.3. An n x n Householder matrix is of the form 

H = I-2vv* (2.1) 

where v is a unit vector. 

Every Householder matrix is a unitary matrix [11, Chapter 5]. Since the 
set of nxn unitary matrices is closed under multiplication, taking the square 
magnitude of entries from a Householder matrix-permutation matrix product 
results in a unistochastic matrix. 

3. Modeling Dynamical Systems 

Consider running a stirring protocol on a compression resistant fluid. Let's 
model this with a measure-preserving dynamical system (D, B, fi, /), 

1. D represents the compression resistant fluid, 

2. i3 is the Borel cx-algebra, 

3. is rescaled Lebesgue measure so that yu(D) = 1, 

4. / : D — ^ D models fluid movement during stirring. 
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The Monte Carlo method we will outline uses nxn unistochastic matrices 
arising from Householder matrix-permutation matrix products. These unis- 
tochastic matrices are close to the permutation matrices when n is large. It is 
reasonable to use the described Monte Carlo distribution for hypothesis test- 
ing when the dynamical system has the following property: For any A,B^B 
where P[A n i?) = 0, if / is perturbed so that 

P{f{x) G A \ X & B) increases (decreases), (3.1) 

then for all Borcl set C C B'^ 

P{f{x) & A \ X & C) decreases (increases) proportionally. (3-2) 

The 1/n row vector and unistochastic matrices arising from Householder 
matrix-permutation matrix products provide Markov shifts that reflect Ulam's 
method with an equal measure partition applied to such dynamical systems. 
This is not saying that all such dynamical systems lead unistochastic ma- 
trices, but that the Householder constructed unistochastic matrices reflect 
these properties. 

If v is a real unit vector and H is the corresponding Householder matrix. 



V2 



H^I- 2v7f. (3.3) 



\Vn/ 

Then squaring the entries of H gives a unistochastic matrix M, 

If Di, Di, . . . , D„ are our equal measure partition sets, and M arose from 
Ulam's method, then the entries of M provide conditional probabilities, 

mij = P{f{x) e I a; e D^)- (3-5) 

Increasing (decreasing) vf leads to nearly proportional decreases (increases) 
in 

rriij = P{f{x) e I X e Di) when i ^ j. (3.6) 

Because of these observations, we propose using unistochastic Ulam matri- 
ces arising from squaring the entries of real Householder matrix-permutation 
matrix products to model weak-mixing stirring protocols of such dynamical 
systems. 
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4. Mixing Hypothesis Test Procedure 

In this section we wiU discuss and outline our procedure for testing a compression- 
resistant fluid stirring protocol. 

There are many techniques to measure a stirring protocol's ability to mix, 
such as decay of correlations [5], Fourier analysis, Artin braid patterns [1, 20], 
chaotic advection [18, 19], and other topological methods [4, 12]. Unfortu- 
nately, it is typical for different protocols to lend themselves to different 
analytical methods. Thus comparing mixing quality between mechanically 
dissimilar stirring protocols is difficult. The primary advantage of the Ulam 
method approximation is that it can be used to evaluate any incompressible 
fluid stirring pattern. This allows one to compare and evaluate the mixing 
of stirring protocols by comparing and evaluating eigenvalues. The main dis- 
advantage is that the method is statistical and does not prove the results. 
Another signiflcant advantage of our method is that it requires only one 
iteration of stirring, in contrast to other techniques that call for iterated 
experiments. 

Since our method only approximates the dynamical system, we make no 
inferences regarding strong-mixing when we conclude that the stirring proto- 
col is weak-mixing. If the protocol is not ergodic, then it is not weak-mixing. 
If the protocol is not weak-mixing, then it is not strong-mixing. 

Our test hypotheses are 

1. Ho : (D, i3,yU, /) is not ergodic (and hence not weak-mixing). 

2. Ha^ : (D,i3,/i, /) is ergodic but not weak-mixing. 

3. '■ /) is weak-mixing (and hence ergodic). 

We partition the fluid into connected, equal volume regions, and use these 
partition sets to generate a new cx-algebra contained in B. If our data strongly 
indicate that the stirring protocol is weak-mixing or ergodic over the gen- 
erated cr-algebra, we will conclude the same about the original dynamical 
system. If the stirring protocol is nonmixing or nonergodic over the gener- 
ated (T-algebra, then the original dynamical system is nonmixing or noner- 
godic. The procedure evaluates stirring over a smaller cr-algebra, thus the 
test is inherently more reliable for detecting if a protocol is nonmixing or 
nonergodic. 

The null hypothesis is that the stirring protocol is nonergodic. It is better 
to reject a protocol that mixes well than to produce poorly mixed product. 
The repercussions of a testing error are as follows: 

1. Type I Error: discard a desirable stirring protocol for a different stirring 
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method 



2. Type II Error: produce a product that is insufficiently mixed 
The stirring protocol's purpose determines the tolerable risks of error and 
the number of partition states. If poorly mixed fluid could result in minor 
consequences or mixing on a small scale is inapt, then the number of par- 
tition regions may be relatively small. If poorly mixed fluid could result in 
severe consequences, then the number of partition regions must be large and 
partition volume small. For example, poorly mixed batter from a kitchen 
could result in unpalatable food; poorly mixed pharmaceuticals with a low 
LD50 could lead to overdose and death. A mixing test for a kitchen could 
use a relatively coarse partition, while a pharmaceutical company would use 
a fine partition. 

If we know an upper bound for the stirring protocol's entropy, call it /i, 
then Froyland's entropy estimate and expected values show that the number 
of states should be greater than [8] . 

Data point movement from one iteration of stirring leads to our empirical 
stochastic matrix, P. The percent of points that start in region i and end in 
region j gives us Pij. We model the entries of P as nonindependent binomial 
random variables, whose probabilities come from the Ulam stochastic matrix. 
Our test statistic is A2 = A2(-P). Some dynamical systems have measure zero 
sets with atypical properties. In an attempt to avoid such difficulties, we 
randomly select data points rather than select points from a grid. 

If the data points are independent, uniform, and randomly distributed 
within each region, then the empirical matrix will converge to a bistochastic 
matrix in a Frobenius norm sense (the proof of this follows from extending a 
standard Monte Carlo argument [16]). We approximate the stirring protocol 

with a one- sided Markov shift, (^{^, . . . , ^) , P^ . The pair will not define a 

Markov shift if 



The 1/n row vector is a stationary distribution for any bistochastic matrix. 

Since Ulam method's partitions our fluid into n equal volume regions, 
it is reasonable to use the l/n vector as the stationary distribution. If we 
do not use equal volume partitions, the stationary distribution will be the 
probability vector corresponding to the rescaled volume of each region. Many 
of the results regarding convergence, expected values, convergence rates, etc. 




n 



1 



(4.1) 
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depend on equal measure partitions, we should use equal volume regions if 
appropriate [17]. 

Mixing Hypothesis Test Procedure: 

1. Set the type II error significance levels for both alternative hypotheses, 
«! and a2- 

2. Set n to be the number of partition regions. 

3. Decide which conditional probability distribution (s) for |A2(-P)| = 1 
and A2(-P) = 1 to use. 

4. Establish the critical values for and ifaj, c\ and C2. The purpose of 
this paper is to propose using Householder matrix-permutation matrix 
products to estimate C\ and C2. 

5. Partition the fluid into n connected equal volume regions, Di, D2, . . . , D„ 

6. Randomly select data points in each partition region. These points 
should be independent and uniformly distributed. 

7. Run the stirring protocol one time. 

8. Use data point movement between regions to 
construct an empirical stochastic matrix, P. 

9. Determine the hypothesis test result. The test statistic is A2(-P); 
compare | A2(-P) — 1 | to Ci; 

compare | A2(-P) | to C2. 
10. Use Froyland's entropy estimate to estimate the dynamical system's 



i=i i=i 

(we define OlogO to equal 0). 
11. If the null hypothesis is rejected in favor of weak-mixing, let the rate 
at which 



entropy. 



n 



n 




(4.2) 




(4.3) 



be our estimate of the rate of mixing. 



5. Constructing the Monte Carlo Matrices 



Let 



nG {3,4,5,6,7,...} 



(5.1) 
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be the number of equal measure states that we partition the measure-preserving 
dynamical system into while using Ulam's method [see [17, Chapter 1] for 
the procedure]. Let 

{Ui,U2, . . . ,Un} (5.2) 

be real independent, identically distributed random variables such that 

Ui^O almost surely, and E{^) , E{u^i), E{\) , E{ut) < oo. (5.3) 

Let 

U2 



U 



\u\- 



(5.4) 



\UnJ 



We may use a unit vector to construct a Householder matrix. Let H = [hij] 
be the Householder matrix corresponding to v. 



H = 1 - 2vv' 



(5.5) 



The entries of H are 

hi T 



1 - 



ru^ if i = j, 



ul+u^ + ...+ul ""i 

^ 2UiUj if i 7^ j. 



(5.6) 



ul+ul+...+ui 



Let Q be a permutation matrix that we want our random unistochastic 
matrix to be proximal to. Set U = QH. Now, let M = [rriij) be the matrix 
defined by 



rriij = Uij. 



(5.7) 



Since i7 is a Householder matrix and Q is a permutation matrix, [/ is a 
unitary matrix. It follows that M is unistochastic. 
Notice that 



rriij 



(1 



u^+ui+...+u 



4 2 2 

iui+ul+...+ul)2^i^j 



if i = j, 
if i ^ j. 



(5.8) 
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Since 7^ almost surely for all z, all entries of M are positive almost 
surely; by the Perron-Frobenius theorem, all but one of M's eigenvalues are of 
magnitude strictly less than one [15, Chapter 8]. It follows that any Markov 
shift with stationary distribution 



will be strong-mixing (for Markov shifts, weak-mixing is equivalent to strong- 
mixing). Our hypothesis test for weak- mixing (ergodic) requires a probability 
distribution over [0, 1] (the unit circle) with significant mass near 1. In the 
next two sections, we will see that the expected value of M's eigenvalues 
converge to one as n goes to infinity. 

How does our dynamical system relate to nl Generally speaking, finer par- 
titions are more apt to detect nonmixing (nonergodicity). If we are confident 
in mixing, we will use a coarse partition to reduce effort; if our confidence in 
mixing is poor, we will use a fine partition. 

By the Birckoff-von Neumann theorem, bistochastic matrices are convex 
combinations of permutation matrices [2, 3]. So our unistochastic matrices 
will tend to be near the 'corners' of the set of bistochastic matrices. For any 
statistic from unistochastic matrices we are interested in, we may use such 
matrices to generate a Monte Carlo empirical probability distribution. 

6. Establishing Critical Values 

In this section, we will outline the procedure we propose for establishing 
critical values for a weak-mixing, ergodic, nonergodic hypothesis test. 
Our test hypotheses are 

1. Ho '■ (D, /) is not ergodic (and hence not weak-mixing). 

2. : (D,B,fi,f) is ergodic but not weak-mixing. 

3. '■ (V>,B,fi,f) is weak-mixing (and hence ergodic). 

After partitioning the space into n equal measure connected subsets, Ulam's 
method approximates the dynamical system with a Markov shift. We will 
approximate the bistochastic matrix defining the Markov shift's measure, P, 
with an empirical stochastic matrix, P. So we approximate 




n 



1 



) 



(5.9) 




(6.1) 



Remark 6.1. The pair 




(6.2) 
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will not define a Markov shift if 

{h---^l)P^{l---^}.)^ (6-3) 
hut if our data points are uniform random variables within each state, then 

E{\\P-P\\f)^Q (6.4) 

as the minimum number of points in a state goes towards infinity. It follows 
that for each eigenvalue 

I Ai(P) — \{P) I—)- in probability V^ (6.5) 

in the Hausdorf topology when our data points are uniform random variables 
within each state and the minimum number of points in a state goes towards 
infinity [17, Chapter 8]. 

Our test statistic is the second largest eigenvalue of P. Let ai, 0^2 G (0, 1) be 
the alpha values for the hypothesis test; let ci, C2 G (0, 1) be the corresponding 
critical values, C2 < 1 — ci, 

P{\ HP) - 1 l> ci : A2(P) = 1) < «i, (6.6) 
P(|A2(P)|<C2:|A2(P)|=l)<a2. (6.7) 

Our goal is to use Householder matrices to estimate Ci and 02- The proba- 
bility distribution used to establish Ci,C2 should reflect properties of a class 
of dynamical systems containing our stirring protocol. 
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Argand diagram of mixing hypothesis test criteria: 

1. If A2(-P) is in the region containing 1 (yellow), fail to reject the 
null hypothesis; conclude that the dynamical system is 
nonergodic. 

2. If A2(-P) is in the outer region away from 1 (green), reject the 
null hypothesis in favor of the first alternative hypothesis; 
conclude that the dynamical system is ergodic, but not 
weak-mixing. 

3. If A2(-P) is in the center region (red), reject the null hypothesis 
in favor of the second alternative hypothesis; conclude that the 
dynamical system is weak- mixing (and hence ergodic). 



Establishing Critical Values for the Test: 

1. Partition D into n equal measure connected subsets. If an upper bound 
of the dynamical system's entropy is known, call the upper bound h, 
set n greater than [8]. 

2. Select a random variable with which to construct unit vectors. 

Let Ml, M2, . . . , Mn be independent, identically distributed, random vari- 
ables, 

M\2- 

3. Set G N so that our empirical probability distributions will be suffi- 
ciently accurate. 

4. Select permutation matrices {Qiji^i near which we want the probabil- 
ity distribution to have significant mass. 

5. Randomly generate Householder matrices. Hi = I — 2viv[. Then 
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square the entries of QiHi to get the matrix Mj. 

6. Use {\2{Mi)}fLi to approximate P{\ X2{P) - I \> k : A2(P) = 1). Use 
the approximation to estimate Ci. 

7. Use {| A2(Mi) to approximate P{\ X2{P) \< k : A2(P) = 1). Use 
the approximation to estimate C2- 

7. Matrix Convergence 

In this section, we will show that Mj from the previous section will converge 
to the permutation matrix Qi as n increases. Since permutation matrix eigen- 
values are on the unit circle, as a random variable, it is likely that the second 
largest eigenvalue from one of our unistochastic matrices will be near mag- 
nitude one. Because of the likely proximity to one, it is reasonable to use a 
probability distribution from such an eigenvalue to establish critical values 
for our weak- mixing, ergodic, nonergodic hypothesis test. 

Our proofs take advantage of the Frobenius norm. After using Jensen's 
inequality to remove the square root from consideration, finding an expected 
value upper bound is similar to finding a second moment. A permutation 
matrix acting on a matrix does not change the magnitude of the entries, 
without loss of generality will prove the results for when Q is the indentity 
matrix and focus on Af s convergence to the identity matrix. 

Proposition 7.1. If M is a matrix constructed in section 2 with Q = I, 
n e {3, 4, 5, . . .}, {ui}2=i are identically distributed, Ui ^ a.s. and 

Eiut),Eiu^),E{^),E{^)<oc, (7.1) 
then E{\\ M — I \\f) as n ^ oo. Moreover, 

m M - I <j^M^)E{u^) (7.2) 
+ J^Ei^)E{ut) (7.3) 
+ '-^E{^mut)r. (7.4) 

(7.5) 

Proof. First, by Jensen's inequality 

E{\\ M - I y) < .Je{\\M-I\\%). (7.6) 
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So it is sufficient to show the second part of the proposition. Let's look at 
the entries of M — 7; by computation we see that: 



-4 



1 



(M-7) 



1-3 



xfl- 



uf+u^+...+u^ 

1 

2 I „2 



(«2+«| + ...+«2) 



if i = j, 
if i ^ j. 



(7.7) 



It follows that 



16 

{ui+ui+...+uir 
32 



{M - 1% = { 



(ul+ul+...+ul)^'^i 



16 



uf if i = J, 



(7i 



If we expand the addends and remove the negative terms, it follows that 
almost surely 



i=l 



16 



i=l 



16 



4 4 

2+ui+...+04^i'"r 



(7.9) 

(7.10) 
(7.11) 



Almost surely, all of the terms in the denominators are positive; if we subtract 
terms from the denominators, we get an upper bound on the fractions. Thus 
almost surely 



M-I ||^< 



16 



1=1 



i=l 



16 



2_i_„2_ 

■2 



16 



.4„.4 



+u: 



2+...+„2_„2_„2)4^i^^iJ■■ 



(7.12) 

(7.13) 
(7.14) 
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Notice that the numerator and denominator in each fraction in the upper 
bound are independent. 

The subtraction of terms in denominators removed some positive terms, 
so the denominators are sums of positive terms. Therefore, we may use the 
harmonic-arithmetic means inequahty, 

II ^ - ^ ll^<16E (^(^ + 4 + ■ ■ • + 4 - ^) (7.15) 

i=l 

n 

1=1 

+ 16y ^^(^ + ^ + ... + 4-^-^)V«t (7-17) 
Now let's take expected values; since the Wj's are independent, 

n 



16 



1=1 



i=l 

Next we use Minkowski's inequality and the fact that the Wj's are indepen- 
dently distributed, 

m\ M - I <j^,E{^)E{u\) (7.21) 

+ J^^E{^)E{<) (7.22) 

+ "^!^E{^){E{<)f- (7-23) 

(7.24) 

Since E{^), E{uf), E(^), and -E(m^) are all finite, it follows that 

E{\\M-I\\l)^0 (7.25) 
as n — >■ oo. Hence by Jensen's inequality, 

E(\\ M - I \\f) ^ (7.26) 

as n — >■ oo. □ 
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Since the second largest eigenvalue gives a test statistic to decide if a 
measure-preserving dynamical system is weak mixing, ergodic, or nonergodic, 
we need a conditional probability distribution of eigenvalues to conduct hy- 
pothesis tests. To statistically test if a measure-preserving dynamical system 
is weak-mixing, we could randomly select permutation matrices 

{Qk}U (7.27) 

and generate {Mk}^^ii with our Householder method, then use the empirical 
probability distribution from 

{\HMu)\}l, (7.28) 

to estabhsh the critical value for the weak-mixing hypothesis test. 

To statistically test if a measure-preserving dynamical system is ergodic, 
we could randomly select permutation matrices 

{Qk ■ the multiphcity of A = 1 is at least two}^]^ (7.29) 

and use Householder matrices to generate {-/V^jtlfeLn then use the empirical 
probability distribution from 

{HM,)}1, (7.30) 
to establish a critical value for the ergodic hypothesis test. 

8. Using Specific Random Veiriables 

In this section, we find more precise upper bounds for specific random vari- 
ables. These upper bounds give better estimates of convergence rate than the 
results in the previous section. The first two proofs in this section start out 
the same way as the first proof in the previous section, then the arguments 
take advantage of the distribution properties. 

Let's find a more precise upper bound when the -Uj's are independent stan- 
dard normal random variables. The proof is similar to the first convergence 
proof, the difference is that we take advantage of the relationship between 
normal random variables and ^^-distributions. 
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Proposition 8.1. // the Ui's in the construction of a unistochastic matrix 
are independent standard normal random variables and n e {11, 12, 13, . . 
then 

m\M-i\\i)< ^^_,,^^}^^:_,,^^_,^ (8.1) 

48ra 
(n-3)(n-5) 

144n 

(n-4)(n-6)(n-8)(n-10) ' 

Proof. Prom the previous proof, we know that almost surely 



+ (n-3Kn-5) f^^-^) 

I 144n(n-l) q\ 

rn-4Vn-6Vn-8Vn-10V 



M - 7 \\l< y , 2^ 2^^i 2 2,4 (8.4) 



1=1 



j=l 

+ V , 2^ 2^ 2 2,4 -^-^l (8.6) 



3 



Since Wj's are independent standard normal random variables, we may 
replace the Wj's with x^-random variables when we take expected values. 



E{\\M-I III) < ^ 16E{j^)E{u^) (8.7) 
1=1 



i=l 



■ (7n-2) 

We use 7j to denote a x^-random variable with j degrees of freedom. If we 
take expected values and remove negative terms, it follows that 

m M - I Kj-,^^,,^^^^ (8.10) 

+ (n-3Kn-5) (^'^l) 

_l 144n(n-l) /o i 9\ 

(n-4)(n-6)(n-8)(n-10) ■ \°--^^) 

□ 
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Now let's consider gamma random variables. A finite sum of independent 
gamma random variables with the same scale parameter is a new gamma 
random variable with the same scale parameter, but the shape parameter 
is the sum of the addend shape parameters. In the next proof, we look at 
independent and indcntically distributed gamma random variables. 

Proposition 8.2. // the Ui 's in the construction of our unistochastic matrix 
are independent T{a,l3) random variables and ^ + 2 < n, then 

Ei\\ M - I <'^^f^f-f^ (8.13) 
^ iy3nLn(^+.) ^g_^4^ 

ni=i[("-l)a-«] ^ ^ 

Proof. Previously we showed that almost surely 

n 

1=1 

n 

+ y" I 2^ 2^^i 2 2,2 (8.17) 

i=l 

^Y. ^^ui^..^-ui-u^^ <^- (8-18) 

Using the Cauchy-Schwarz inequality and the fact that < tij almost surely 
for all i, we get 

n 

\\M-I \\l< y r^-^ ^u^ (8.19) 

1=1 

n 

+ y" f -u (8.20) 

/ ^ (ui+U2+—+Un-UiP « V / 

i=l 

+ Y < ^ ^ (8-21) 

If we take expected values, and take advantage of the independent and iden- 
tically distributed lij's. 
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EiW M - I 11^) <^( (.,,.„,;!t„_.,)s )^(^^D (8.22) 

' 16nf 

^ {Ul+U2 + ---+Un-Ui)* 

' 16nS(n-l) 

^ {Ul+U2 + ---+U„-Ui-Uj)^ 



+ E( ,^..uJ'rL-u,M ut) (8.23) 
+ E L..Ji!^l-u.^s )Eiv^)Ein^)- (8.24) 



Since the Wj's are gamma random variables and ^ + 2 < n, wc may replace the 
denominator sums with gamma random variables. Let 7j denote a gamma 
random varible with parameters ja, and /3. 

£(||M-/|||)< g;n;;;^^ (8.25) 

□ 

Next we will consider bistochastic matrices arising from random House- 
holder matrices where 



ui — U2 — ■ ■ ■ — Un, and Ui (8.28) 
almost surely. It follows that almost surely 



For the previous matrix probability distributions described, we use random 
numbers to build matrices, then use the matrices to construct an empirical 
approximation to the probability distribution of the statistic under consid- 
eration. When Ui = U2 = ■ ■ ■ = Un, Ui ^ almost surely, the aspects of the 
matrix only depends on n, we know the precise distribution of all matrix 
statistics. If reasonable for the dynamical system under consideration, we 
could use M to estabhsh critical values for hypothesis testing. So let's look 
at some aspects of such M's, after a lemma. 



,M=('^)/ + 



(8.29) 
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Lemma 8.3. // Dn and Sn are n x n symmetric matrices such that 





a 




P . 


■ a 








OL 


P ■ 


■ P 








P 


a . 


■ P 


and Sn 




/3 




a ■ 


. a 





^ Dn-1 



where /3 is the n — 1 vector with /3 for all entries, then 

dct(AO = {a- /3T'\a + (n - and 
det{Sn) = {a- py-^p. 



(8.30) 



(8.31) 
(8.32) 



Proof by induction:. Assume true for n. Using the fact that interchanging 
any two rows or any two columns of a real matrix changes the sign of the 
determinant, we see that 



det{Dn+i) =Q;det(D„) — ^ndet(5'„) and 
det(,S„+i) =^det(L'„) -^ndet(,S„). 

Using the induction hypothesis, these equations become 



(8.33) 
(8.34) 



det(D„+i) =a{{a- /?)"-'(« + (n - l)/3) ] - l3n{ {a - l3)''-^l3 and 



.35) 



det(>S„+i) \^(a - p)''-\a + (n - - pn \^(a - PT'^Pj . (8.36) 

Factoring out the (a — /3) terms gives us the results. □ 
Proposition 8.4. If M is the n x n matrix matrix 



then det(M) = (^)" \ trace{M) 
of M is the diagonal matrix with entries (1, 



(n-2)2 



(8.37) 



and the Jordan canonical form 



n— 4 



n— 4 



)• 
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Proof. The trace of M follows from the definition. The determinant follows 
from the previous lemma by setting a = ^"'~2^ , and /? = 

Now, the matrix is a symmetric real matrix; hence it is diagonalizable. The 

eigenvalues follow from the previous lemma by setting a = ^"~2^^ A, and 

(3 = ^2 to get the characteristic polynomial of M. □ 

Since M is diagonalizable, if the Markov shift {{l/n,l/n, . . . ,1/n), M) 
approximates a measure-preserving dynamical system that is mixing, the 
estimate of mixing rate is the rate at which 

(^)^^OasiV^oo, (8.38) 

instead of the estimate given in [17, Chapter 4], (^7^)^""^^ — )• as 
N ^ 00. 



9. Two Region Partitions 



There are few instances of interest where one would use our method with two 
partition regions. We look at this special case as an example to help develop 
understanding. A potential application is equal ratio mixing of items with 
minimal consequences of poor mixing, such as combining blends of coffee. 
When combining two equal volumes of coffee, poor mixing would result is 
inconsistent taste. Only the most serious baristas would say that inconsistent 
cup-of-Joe flavor is worse than poorly mixing pharmaceuticals. 
Say that our unit vector is 



^^1 

Since the vector has norm one, we may write its Householder matrix as 
H 



(9.1) 



1 - 2vl 



±2^1 VI 



(9.2) 



I 2vf - 1 

There are two possible doubly-stochastic matrices arising from our method. 



and 



:i-2vir 



(9.3) 
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Whose characteristic polynomials are 

(1 - X){8v^ - 8vl + 1 - A), and (1 - A)(-8ff + Sv^ - 1 - A). (9.4) 

The second largest eigenvalue depends on Vi. Let's graph the relationship. 

-A2 

h 




Graphs of the relationship between vi and A2 when n = 2 

When n = 2, the relationship between vi and second largest eigenvalue 
defines a function from [—1,1] to [—1,1]. In this situation, it is feasible to 
compute 

P(|A2 - 1| > A; : A2 = 1) and P(|A2| > k : IA2I = 1). (9.5) 
If f 1 is a beta random variable with parameters a and /3, then 



^(A2 



± 



8a(a+l)(a+2)(a+3) 



(a+/3)(a+/3+l)(a+/3+2)(a+/3+3) 



8a(a+l) 



± 1, 



(9.6) 
(9.7) 



^ (a+^)(a+/3+l) 

where the sign of addends depends on the permutation matrix used. 
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